clear all
set more off
set scheme s2color

// Density of network around regional borders

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_data.dta", clear

ren com_res com
collapse (firstnm) dep_res (mean) dist_border (sum) ind, by(com)

****************** Costmetics ******************
** Remove Corsica
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com, replace


****************** Get dist border ******************
merge 1:1 com using "$input\com_arcgis.dta"
drop if _merge != 3
drop _merge
** Strict on natural features relative to rest. 
keep if dist_nat_border < 0 | dist_nat_border > 2
keep if dist_river < 0 | dist_river > 10

****************** Gen length network ******************
merge 1:1 com using "$input\road_dep_intersect.dta"
drop if _merge != 3
drop _merge

****************** Get neighbor region for all com ******************
tostring com, replace
replace com = "0"+com if strlen(com) == 4
merge 1:1 com using "$input\loc_i_raw.dta" 
drop if _merge == 2
drop _merge
destring com dep_res, replace

********* Other side ****************
geonear com x_com y_com using "$input\loc_j_raw.dta" , neighbors(dep x_dep y_dep) nearcount(2)
gen other_dep = nid1 if nid1 != dep_res
replace other_dep = nid2 if nid1 == dep_res

egen home_lab_mkt = sum(ind), by(dep_res)
egen other_lab_mkt = sum(ind), by(other_dep)


****************** Gen vars ******************

replace dist_dep = -dist_dep if other_lab_mkt > home_lab_mkt

gen r_dist_dep = round(dist_dep + 1.01, 2) if other_lab_mkt < home_lab_mkt
replace r_dist_dep = round(dist_dep - 1.01, 2) if other_lab_mkt > home_lab_mkt

keep if abs(r_dist_dep) < 16

replace r_dist_dep = r_dist_dep + 16
tab r_dist_dep

gen length_pkm2 = length/surface
gen length_pc = length/ind

reghdfe length i.r_dist_dep x_com y_com, absorb(dep_res other_dep)
reghdfe length_pkm2 i.r_dist_dep x_com y_com, absorb(dep_res other_dep)
est store reg_1

coefplot (reg_1,  color(black) ciopts(recast(. rcap) color(black black) ///
lwidth(*0.75 *1.25)) ), base ///
vertical omitted keep(*.r_dist_dep*)  yline(0) graphregion(color(white)) bgcolor(white) ///
msize(small) levels(95 90) xline(7.5) ///
xlabel(1 "-14{&le}.<-12km" 2 "-12{&le}.<-10km" 3 "-10{&le}.<-8km" 4 "-8{&le}.<-6km" 5 "-6{&le}.<-4km" 6 "-4{&le}.<-2km" 7 "-2{&le}.{&ge}0km" 8 "0{&le}.<2km" 9 "2{&le}.<4km" 10 "4{&le}.<6km" 11 "6{&le}.<8km" 12 "8{&le}.<10km" 13 "10{&le}.<12km" 14 "12{&le}.<14km",  angle(45)) ///
ytitle("Road density (km/km2)") ///
text(0.035 5.3 "Departmental border", size(Large) box fcolor(none) lc(green)) ///
text(-0.0575 5.7 "{bf:Smaller market}", size(Large) ) ///
text(-0.0575 9.2 "{bf:Larger market}", size(Large)  ) 
graph export "$output\network_gap_formal.pdf", replace

